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The  differential  system 


Consider  the  matrix  differential  equation 


Y(t)  =  Y(tyTF(Y(t),Y(tyT) 
¥(())  =  70  e  GL(n) 


F  is  a  continuous  matrix  function,  globally 
Lipschitz  on  a  subdomain  of  GL(n) 

the  solution  Y(t )  exists  and  is  unique  in  a 
neighborhood  ]-t  t[  of  the  origin  0 
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The  structure  of  GL(n) 


♦>  Two  maximal  connected  and  disjoint  open  subsets 
comprising  GL(n) 
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♦> 


Theoretical  results 


The  existence  of  the  solution  Y(t)  for  all  t  is  not 
guaranteed  a  priori  and  the  presence  of  a  finite 
escape  time  behavior  is  not  precluded. 

The  value  of  the  escape  point  depends  on  the 
function  F 

>  If  the  escape  point  r  is  finite  then  Y(t)  approaches 
singular  matrix  as  t  — >  z 

y  if  t  <  oo  then  Y(t)  exists  for  all  t  >  0 
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Theoretical  results 


Example:  F  constant  function  with  trace{F)  =  0 


Y  =  Y~t 


0 

-1 


j_  VT+7 

a/2  VW 
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Theoretical  results 


❖  Relationship  between  the  singular  values  of  the 
solution  Y(t),  the  initial  condition  Y(0 )  and  the 
symmetric  matrix  function: 

t 

E(t )  =  j[Fr(FCs-),rr(,y))+ F(Y(s) ,  Y~r  (.s))]  ds 

0  _ 


Systems  with  structure 


If  the  matrix  function  F  maps  all  matrices  into  the 
Lie  algebra  of  skew-symmetric  matrices 

Y(t)  belongs  to  the  orthogonal  manifold 
(whenever  Y(0 )  is  orthogonal) 

If  diag(F)  =  0  for  all  nonsingular  matrices 
diag(Y(t)TY(t))  =  diag(Y(0)TY(0)) 
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Examples 


Control  Theory 

>  Optimal  system  assignment  via  Output  Feedback 
Control 

>  Balanced  Matrix  Factorizations 

>  Balanced  realizations  (Isodynamical  flows) 

Multivariate  Data  Analysis 

>  Weighted  Obliqne  Procrustes  problem 

Inverse  Eigenvalue  Problem 

>  Pole  placement  or  eigenvalue  assignment  problem 
via  output  feedback 

>  Prescribed  Entries  Inverse  Eigenvalue  Problem 
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Examples  in  Control  Theory 


Output  Feedback  Control  of  linear  system 


>  Consider  the  linear  dynamical  system  defined  by  the 
triple  (A,B,C)e  p«  x  «  X  Pw  x  w  X  P/7  x  w 

x(t )  =  Ax  ( t )  +  Bu  (t) 
y(t)  =  Cx(t ) 


>  The  process  of  “feeding  back V  the  output  or  the  state 
variables  in  a  dynamical  system  configuration  through 
the  input  channels 


>  Output  Feedback:  u(t )  is  replaced  by  u(t)=Ky(t)+v(t) 


£ej>mxP  feedbaCk  gain  matrix 
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Examples  in  Control  Theory 


Output  Feedback  Control  of  linear  system 

>  The  feedback  system  is 

x(t)  =  (A  +  BKC  )jc(0  +  Bv  (t) 
y(  0  =  Cx(t) 

Optimal  system  assignment 

>  Given  a  target  system  described  by  the  triple 

(F,G,H)eFnXnx¥nXmx¥pXn  find  an  optimal  feedback 
transformation  of  (A,B,C)  which  results  the  best 
approximation  of  ( F,G,H ). 
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Examples  in  Control  Theory 


♦♦♦  The  set  GL(n)xPmXP  of  feedback  transformation 
is  a  Lie  group  under  the  operation 

(T„K,)o(T2,K2)  =  (T,T2,K,+K2) 

♦♦♦We  can  consider  action  on  the  output  feedback 
group  and  orbits,  particularly: 

<I>(A,B,C)  =  {(T(A  +  BKQT~\ TB, CT~ 1  \(T,K)<=  GL(n ) x P mxp 


♦♦♦  The  distance  function 
®  =||  T(A  +  BKC)T~l  -  F 


+  TB-G 


CT~l  -  H 
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Examples  in  Control  Theory 


♦>  The  gradient  flow  of  this  distance  function  with 
respect  to  a  specific  Riemannian  metric  on 
<t>(A,B,C)  can  be  written  as: 

T  =  T~t  f(T,T~r  ,K) 

K  =  -BtTt(T(A  +  BKC)T~l  -  F)T~tCt 


Nicoletta  Del  Buono 


Examples  in  Control  Theory 


Balanced  matrix  factorizations 

>  General  matrix  factorization  problem: 

Given  a  matrix  HePkxl find  two  XePkXn  and  YePnXl 
such  that  H=XY 

>  balanced  factorization  XTX=YYT 

>  diagonal  balanced  factorization  X TX=YYT=D 


Balanced  and  diagonal  balanced  factorization  can  be 
characterized  as  critical  points  of  cost  functions 
defined  on  the  orbit 


O  (X,  7)  =  {(. XT ■' ,  TY)  e  P  kxn  x  P  ”xZ 


T  e  GZ(n)} 
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Examples  in  Control  Theory 


❖  The  cost  functions  are  respectively: 


<t>  :0(X,y)  — »  P 
0>„:O(X,y)->P 


<X>(Xr\TY)  =||  XT 


-1  ||2 


O  TXT1 ,  TY )  =  tr(NT~l  X1  XT'  +  NTYY’  T‘ ) 


+  ||  TY 

~T  t rT 


T-l 


'TrrT 


❖  Applying  a  gradient  flow  techniques  differential 
systems  on  GL(n)  can  be  constructed: 

iced  • 

-*~T  =  T~t (XT X(TtT)~x -TtTYYt)  T(0)  =  T0 

T  =  T-t(XtXT-1NT-t  -TtNTYYt)  T(0)  =  7;^^7 

balanced 
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Examples  in  Control  Theory 


Balanced  realizations  in  linear  system  theory 

>  Consider  the  linear  dynamical  system  defined  by  the  triple 

(A  AQ  eP  nxnxpnxmxppxn 

x(t)  =  Ax  ( t )  +  Bu  ( t ) 

y(t)  =  Cx(t ) 

00  00 

>  Gramians:  Wc=  J  eMBBTeAftdt  W0  =  j  eATtCTCeA,dt 

0  0 

>  (A,B,C)  is  a  balanced  realization  if  Wc=  WQ 

>  (A,B,C)  is  a  diagonal  balanced  realization  if  WC=W0=D 
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Examples  in  Control  Theory 


Any  T e  GL(n )  changes  a  realization  by 
( A,B,C)^(TAT TB,CT~!) 
and  the  Gram  i  an  s  via 

WC-^T  Wc  T~l  l¥0^T-TW0T-! 


Balanced  and  diagonal  balanced  realizations  have 
been  proved  to  be  critical  points  of  costs  functions 
defined  on  the  orbit 


O  (A,  B,C)  =  {(! TA  r 1 ,  TB,  Crl  )eP“xP“xP ,fX" 


T  e  GL(n )} 
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Examples  in  Control  Theory 


❖  The  cost  functions  are  respectively: 


❖  All  balancing  transformation  T e  GL(n )  for  a  given 
asymptotically  stable  system  (A,B,C)  can  be 
obtained  solving  the  gradient  flow 

balanced  J  =  J-T  (Wq(TT  T)~l  -  TT  TWc)  T  (0)  =  TQ 

T  =  T  ~T  {W 0T NT  ~T  -  Tr  NTWC)  T{<d)  =  Ta  diagonal 

balanced 


$:0(i,S,C)^P  0(T)  =  tr(TWcT~l  +  ) 

<Z>N  :  O  ( A,  B,  C)  ->  P  O  ,,(T)  =  tr(NTWcT +  NT~tW0T~' 
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Examples  in  Multivariate  Data 

Analysis 


Weighted  oblique  Procrustes  problem  (WObPP) 

>  Manifold  of  the  oblique  rotation  matrices 

OB(n)  =  {JeP“|  det(X)  *  0  ,diag(XT  X)  =  1} 

Given  A,B,  C  fixed  matrices  with  conformal  dimensions 

>  Minimize  ||  AXC-  B  ||  subject  to  XeOB(u) 

>  Problem  in  factor  analysis  known  as  a  “rotation  to 
factor-structure  matrix” 

>  Minimize  ||  AX  TC-  B  ||  subject  to  XeOB(n) 

>  Problem  of  finding  an  approximation  to  a  “ factor - 
pattenrm&tnx 
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Examples  in  Multivariate  Data 

Analysis 


The  solution  of  the  WObPP  problem  can  be 
obtained  solving  a  descent  matrix  ODE: 


dX 

dt 


=  -^(,)(V)  =  -^‘ro#(^rV) 


being  V  the  gradient  of  the  function  to  be 
minimize  with  respect  to  the  chosen  metric 

(N.  Trendafilov  FGCS  2003) 
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Examples  in  Inverse  Eigenvalue 
Problem  and  control  theory 


Pole  placement  or  eigenvalue  assignment  via 
output  feedback: 

>  Given  a  linear  system  described  by  the  triple  (A.B.C) 
and  a  self-conjugate  set  of  complex  points  I/.,/., ...  X  } 

>  find  a  feedback  gain  matrix  K  such  that  A+BKC  has 
eigenvalues  /.,■ 

Denoted  by  A  a  fixed  matrix  with  eigenvalues 
the  pole  placement  task  is  equivalent  to  find  a 

matrix  TeGL(n)  and  ATeP  mXP  minimizing  the 
distance  pi  y  A_ T^A  +  BKC)T~X  \ \ 
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Examples  in  Inverse  Eigenvalue 
Problem  and  control  theory 


❖  Using  a  gradient  flow  techniques  the  solution  can 
be  obtained  solving 


T  =  T~t[{A  +  BKC  )T ,Tt (h  -  (A  +  BKC  ))T~T] 
K  =  -BtTt(T(A  +  BKC  )T~l  -  F)T~tCt 
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Examples  in  Inverse  Eigenvalue 

Problem 


❖  Matrix  completion  with  prescribed  eigenvalues 

❖  PEIEP  (prescribed  entries  inverse  eigenvalue 
problem) : 


Given 

>  A  ={(ivJy)  |  v  =1  m  pairs  of  integers  1  <iv  <jv<n 

>a={al,...,am}  c  P 

>  { Ip ln }  a  X  closed  under  conjugation 

Find  a  matrix  A  €  P"x"  such  that  a(X)={  /., . } 

and  x..  =a  v  =  l 

l vJv  ” 
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Examples  in  Inverse  Eigenvalue 

Problem 


►  Let  A  a  matrix  with  eigenvalues  Xt  and  denoting 

M(A)={VAV~l  |  VeGL(n )  } 

the  orbit  of  matrices  isospectral  to  A  under  the 
action  group  of  GL(n)  and 


S(A,«)={X=fc]s  P 


nxn 


x  =  a  v  =  L... 

/  7  v  v 


vJv 


mj 


Solving  the  PEIEP  is  to  find  intersection  of  the 
two  geometric  entities  M(yf )  and  E(A,«) 
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Examples  in  Inverse  Eigenvalue 

Problem 


❖  Minimize  for  each  given  Xs  M(A  )  the  distance 
between  X  and  X(A  ,a) 

min  -<  VAV~l  -  P(VAV~l),VAV~x  -P(VKV~l)  > 

FeM(A)2  | 

Projection  on  I(A,a) 

❖  Using  a  descent  flow  approach  we  get 

—  =  k(VKV~1)V~t  with  k(X)  =  [X\X-P(X)\ 
dt 

(M.T.  Chu  et  al.  FGCS  2003) 
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Numerical  Approximation: 
substituting  approach 


*>  Consider  our  system: 


Y(t)  =  Y  (tyT  F(Y  (0,  Y  (tyT) 
7(0)  =  70  e  GL(n) 


♦>  Setting  Z=Y  ~T  from  YTZ=I  we  get 


J 

Y  =  ZF(Y,Z)  =  H(Y,Z), 

Y(0)  =  Y„ 

-s 

• 

Z  =  -ZFt{Y,Z)ZtZ  =  -ZHt{Y,Z)Z, 

Z(0) =  y;t 

Nicoletta  Del  Buono 


♦> 


Substituting  Approach 


Advantages: 

>  No  direct  use  of  the  inverse  of  Y(t )  (computational 
advantages) 

♦♦♦Drawbacks: 

>  Solution  of  a  new  matrix  ODE  with  doable  dimension 
with  respect  to  the  original  system; 

>  High  stiffness  (when  Y(t)  tends  to  a  singular  matrix  or 
the  Lipschitz  constant  of  H  is  large); 

>  The  presence  of  an  additional  structure  of  the  solution 

matrix  Y(t)  is  not  considered  need  of  ad  hoc 

numerical  scheme 
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Solution  via  Riccati  equation 


When  the  matrix  function  F  does  not  depend 
explicitely  on  Y~T,  i.e.: 


Y(t)  =  Y(tyr  F(Y(t)) 
7(0 1  =  70  e  GL(n) 


►  It  could  be  convenient  work  with  the 
equation 


Y  {t)Y  (t)  =  F  (Y  (t)) 
7(0)  =  70  e  GL(n) 
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Solution  via  Riccati  equation 


*>  Applying  the  second  order  Gauss  Legendre 
method,  we  get: 

Y+Y 


C.r_ , +y;ym -y;j„-y;y„  -2afi 

V 


A 


n+ 1 


=  0 


J 


*>  The  previous  equation  can  be  iteratively  solved 
starting  from  an  initial  approximation  Yn(0) 

(avoiding  the  nonlinearity  of  F) 


n+ 1 


f 


Y 1  Y  +7  7  -Y  Y 

1 71+1 1 72+1  ^  1  n  1 22+1  1  n+l1  n 


YX-2hF 


Yn 


(01\ 


n+ L 


V. 


V 


=  0 


J 
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Solution  via  Riccati  Equation 


The  latter  equation  is  the  prototype  of  an  Algebraic 
Riccati  equation,  in  fact  setting 


A  =  Y„  and  C  =  Y  Y  +  2hF 


Y  +Y 

n  r 


(0)\ 


72+1 


V 


J 


we  get 


R(X)  =  XrX  +  AT  X  -  XT  A  +  C  =  0 
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Solution  via  Algebraic  Riccati 

equation 


Numerical  methods  to  solve  Algebraic  Riccati 
equation  are  based  on  fixed  point  or  Newton 
iteration: 


>  Picard  iteration: 


>  Newton  method: 


n*n 


inxn 


>R:  P 

>  its  Frechet  derivatitive  is:  RX(H)  =  HT (X  -  A)+(X + A)T H 

>  the  Newton  iteration  starts  from  X0  and  solves 
7?(X)=0  via  Xk+1=  Xk+  Dk  being  Dk  the  solution  of 


Sylvester  equation 
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Solution  via  Riccati  equation 


♦♦♦  Solving  Riccati  equation  implies  the  numerical 
treatment  of  the  Sylvester  equation 

AX+XTB  =  X 


with  A,  B,  X  given  n  xn  matrices 


Existence:  there  exists  a  solution  X of  the 
Sylvester  equation  iff 


X  A 
B  O 


A 

O 


are  equivalent 
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Solution  via  Riccati  equation 


♦♦♦To  obtain  conditions  for  uniqueness  of  solution 
and  for  constructing  it,  we  reformulate  the 
Sylvester  equation  as  a  n2n2  linear  system: 


(7  ®  A)vec(^)  +  (Bt  ®  I)vec(Xr)  =  vec(X) 


vec(XT)  =  P(n,n)vec(X) 
P(n,n)  =  ±±Eij®E!J 

i= 1  7=1 


[(7®  A)  +  (Bt  ®  I)P(n,n)]  vec(X)  =  vec(X) 

V - v - J 

M 
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Solution  via  Riccati  equation 


A  +  eyb[ 

ebT 

^2U\ 

e  bT 

M  - 

ebT 

^\U2 

• 

• 

A  +  e2bl 

•  • 

•  • 

e  bT 

^nu 2 

• 

• 

ebT 

L  ^1  Un 

etf  - 

A  +  enbTn 

being  bt  the  columns  of  the  matrix  B 

Uniqueness:  there  exists  a  unique  solution  X  of  the 

Sylvester  equation  AX+X  TB=X  if  the  matrix  M  is 
non-singular  (rank  ( M)=n 2) 
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Solution  via  Riccati  equation 


Considering  the  linear  equation  derived  from: 

>  Picard  iteration:  A =AT  and  B  =A  M is  singular 

>  Newton  iteration:  A  =Xk+AT  and  B  =Xk-A  =>  M  is 
non-singular  =>  unique  solution  ! 

Newton  method  converges  in  a  reasonable  number 
of  iterations 

Numerical  solution  of  Sylvester  equation  : 

>  Direct  methods  (QR„  Gaussian  Elimination); 

>  Iterative  algorithms; 

>  Generalize  Conjugate  Residual  method. 
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Singular  Value  Decomposition 


To  avoid  the  inverse  matrix  computations  and  to 
control  the  singularities  of  the  matrix  solution  Y(t)  we 
can  adopt  a  continuous  Singular  Value  Decomposition 
approach 

The  continuous  SVD  of  Y  (t)  is  a  continuous 
factorization 

m=m  mvr(t) 

>  U(t),V(t )  ortliogonal  matrices  I  U1 U=I„  and 

y  Z(t)  diagonal  matrix  with  diagonal  elements  the 
singular  values  oi  (t)  of  Y(t) 

The  motion  of  Y(t)  is  now  described  by  the  variables  U(t), 
V(t)  giving  more  information  on  the  flow 
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Singular  Value  Decomposition 


❖  Suppose  that  the  solution  Y(t)  possesses  dinstinct 
and  nonzero  singular  values  <Jj(t),  for  i=l,...,  n 

and  t  in  [0,  x)  then  there  exists  a  continuous  SVD 
of  Y(t)  and  the  factors  U(t),  E(t),  V(t)  of  such  a 
decomposition  satisfy  the  following  ODEs: 

i  =  X-1VTF(Y,Y-T)F  -HY  +  YK,  £(0)  =  £0 
U  =  UH,  U( 0)  =  U0 
V  =  VK,  V(0)  =  V0 
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Singular  Value  Decomposition 


♦♦♦  The  differential  equations  for  the  singular  values  are 
=  —tyTF(YJ-T)v\,  i  =  1,-  ,n 

(7  ■ 


♦>  The  elements  of  the  skew-symmetric  matrices  H,  K 
are 


H„ 

1 

- T~  ^  / 

(vtFV  1  +  <r 2  (v  T FV  ).  1 

V 

<r?)L  J 

V  Ay  i  \  Aji  J 

K„ 

=  '  1 

\(v  T  FV  ) 

,  +  (v  1  FV  )  1 

V 

K2-^,2)1 

LV  A 

J  \  Aji  J 
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Singular  Value  Decomposition 


Numerical  solution  of: 

>  a  diagonal  system  in  ai  (information  on  the  conditioning 
of  the  matrix  solution  Y{t)) 

>  two  linear  systems  in  H-  Kl} 

>  two  orthogonal  systems  in  U  and  V 

>our  aim  is  to  preserve  the  non-singular  behavior  of 
the  numerical  solution  ->  explicit  integration  of  the 
systems  in  V  and  V  (orthogonality  preserved  up  to  the 
order  of  the  method) 

Drawback  >=♦  distinct  singular  values 

>  Block  Continuous  SVD 
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Rectangular  case 


Some  of  the  previous  results  can  be  extended  to 
differential  problems  on  the  manifold 

GL(m,  n)  ={YePmXn  |  rank(Y)  =  n},  n  <  m 

Differential  systems  on  GL(m,ri)  have  the 
following  form: 


Y  =  G{Y),  Y (0)  -Y0e  GL{n, p) 


with  G  belonging  to  the  tangent  space  of  GL(m,n ) 


G(Y)  =  r(rT¥) 


+  i,-y(y’yYyt 
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Rectangular  Case:  numerical 

treatment 


Continuous  SVD  (economy) 

777  x/7 

r-^ 

II 

/77  x  /7  matrix 

nxn  matrix 

u,Tut=in 

diag(a a n) 

VTV=VVT_=In 

Differentiating  we  obtain  the  differential  systems 
satisfied  by  the  three  factors: 
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Rectangular  Case:  numerical 

treatment 
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Rectangular  Case:  numerical 

treatment 


Substituting  approach: 


Y  =  y{jtyY  Fx(Y)  +  [/„  -  y{ytY )  Fr]  F2(Y) 


Setting  Z=(Y  TY  )_1  we  obtain 

Y  =  YZF,(Y)  +  [/  -  YZY  r]/s(}') 
Z  =  -Z[F1(Y)  +  FIT(Y)\Z 
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♦>  ♦> 


Numerical  Illustrations 


First  example: 


Y  =  Y~t 


0 

S 

2 


S 

2 

0 


With  solution  existing  in  (-1/3, 1/3) 


J_  Vl  +  St 

42  4i-st 


-4\+st 

4\-st 


We  solve  the  problem  with  <5  =1/2 
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Global  Error 


Numerical  Illustrations 


♦>  Behaviour  of  the  global  error  on  [0  2) 
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♦>  ♦>  ^  .  ♦> 


Numerical  Illustrations 


Second  example 


-Y~t 


sin(f)  cos(t)  cos(f) 
-/sin(/)  t 


1  0 
0  1 


with  solution 


cos(t)  t 
0  1 


periodically  singular  (for  each  xk=  k  n! 2) 
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Global  Error 


Numerical  Illustrations 


❖  Semilog  plot  of  the  global  error  on  {n! A, Till) 
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♦>  ♦>  ♦> 


Conclusions 


We  have  considered  a  particular  ODEs  on  GL(n) 
often  occurring  in  applications 

Several  problems  modeled  by  such  ODEs 

Different  numerical  approaches  avoiding  the 
direct  use  of  matrix  inversion  and  detection  of 
singular  behavior 

Future  works: 

>  Improving  the  validation  of  the  proposed  approaches  by 
tackling  numerical  tests  on  real  examples 
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